uniform sampler2D texUnit;
uniform mat4 direction;
uniform float dtime;
uniform int isOpacity;

float torus(in vec4 inVec)
{
	float x2, y2, z2, f;
	vec4 testVec = vec4(inVec.x, inVec.y, inVec.z, 1.0);

	x2 = testVec.x*testVec.x;
	y2 = testVec.y*testVec.y;
	z2 = testVec.z*testVec.z;
	f = 0.09 - x2 - y2 - z2 - 0.25 + sqrt(x2+y2);
	return f;
} 


bool func(in vec4 inVec, inout float fValue)
{
	fValue = torus(inVec);
	if (fValue >= 0.0)
	{
		return true;
	}
	else 
	{
		return false;
	}

}

bool find_point(inout vec4 vecBegin, in vec4 vecEnd, in float fIntervalLength, in int nReativePrecision, in float fFuncPrecision)
{
	int nIt, nIterations, nSubIt, nSubIterations;
	vec4 vecDiff, vecPoint;
	vec4 vecIntervalBegin, vecIntervalEnd, vecIntervalMiddle;
	bool bMiddleTest, bReturn, bBegin, bEnd;
	bReturn = false;
	float fDistance, fBegin, fEnd, fDummy;

	fDistance = length(vecEnd - vecBegin);
	nIterations = int(fDistance/fIntervalLength);
	nSubIterations = nReativePrecision;
	bBegin = func(vecBegin, fBegin);
	if (bBegin){
		nIterations = 0;
	}

	vecDiff = normalize(vecEnd - vecBegin);
	vecDiff *= fIntervalLength;
	vecPoint = vecBegin;
	fDistance = fIntervalLength;
	fEnd = fBegin;
	for (nIt = 0; nIt < nIterations; nIt++)
	{
		fBegin = fEnd;
		vecPoint += vecDiff;
		bEnd = func(vecPoint, fEnd);
		if (bEnd)
		{
			//r_intersection found
			vecIntervalBegin = vecPoint - vecDiff;
			vecIntervalEnd = vecPoint;
			vecIntervalMiddle = (vecIntervalBegin+vecIntervalEnd)*0.5; 
			nIt = nIterations;
			bReturn = true;
		}
	}

//chord method
	for (nSubIt = 0; nSubIt < nSubIterations; nSubIt++)	
	{
		vecIntervalMiddle = vecIntervalBegin - (vecIntervalEnd - vecIntervalBegin)*fBegin/(fEnd-fBegin);
		bMiddleTest = func(vecIntervalMiddle, fDummy);
		if (bMiddleTest)
		{
			vecIntervalEnd = vecIntervalMiddle;
			fEnd = fDummy;
		}
		else
		{
			vecIntervalBegin = vecIntervalMiddle;
			fBegin = fDummy;
		}
		if (abs(fDummy) < fFuncPrecision)
		{
			nSubIt = nSubIterations;
			break;
		}
	}
	
	vecBegin = vecIntervalMiddle;
	return bReturn;
}

bool find_point_proxy(inout vec4 vecBegin, in vec4 vecEnd, in float fIntervalLength, in int nReativePrecision, in float fFuncPrecision, out vec4 vecIsec)
{

	bool b;
	b = find_point(vecBegin, vecEnd, fIntervalLength, nReativePrecision, fFuncPrecision);
	vecIsec = vecBegin;
	return b;
}

vec3 find_normal(in vec4 vecPoint, in float fPrecision)
{
     vec4 dx, dy, dz;
     dx = vec4(vecPoint.x + fPrecision, vecPoint.y, vecPoint.z, 1.0);
     dy = vec4(vecPoint.x, vecPoint.y + fPrecision, vecPoint.z, 1.0);
     dz = vec4(vecPoint.x, vecPoint.y, vecPoint.z + fPrecision, 1.0);
     float testX, testDX, testDY, testDZ;
     func(vecPoint, testX);
     func(dx, testDX);
     func(dy, testDY);
     func(dz, testDZ);
     return normalize(vec3(testX-testDX, testX-testDY, testX-testDZ));
}

void main(void)
{
	vec2 texCoord = gl_TexCoord[0].xy;
	vec4 c, bbox_min, bbox_max, pos, pend, pend_tr, pos_tr, isec_tr, z;
	vec3 norm_tr;
	bool b1, b2, b;
	int iter;
	float fStep, fPrecision;
	float fDepth, fDepthLength;

	//bounding box by default
	if (isOpacity == 0)
	{
	bbox_min = vec4(-1.0,-1.0,-1.0,1.0);
	bbox_max = vec4(1.0,1.0,1.0,1.0);
	}
	else
	{
	bbox_min = vec4(-1.0,-1.0,1.0,1.0);
	bbox_max = vec4(1.0,1.0,-1.0,1.0);
	}
	//obtaining position
//	c = texture2D(texUnit, texCoord);
	z = vec4(0.0,0.0,0.0,0.0);
	b = false;
	c.x = texCoord.x;
	c.y = texCoord.y;
	c.z = 0.0;
	c.w = 1.0;
	pos = vec4(bbox_min.x + (bbox_max.x - bbox_min.x)*c.x, bbox_min.y + (bbox_max.y - bbox_min.y)*c.y, bbox_min.z + (bbox_max.z - bbox_min.z)*c.z, c.w);
	pend = pos + vec4(0.0, 0.0, (bbox_max.z - bbox_min.z), 1.0);
	pos_tr = vec4(direction * pos);
	pend_tr = vec4(direction * pend);
	
	fStep = length(pend_tr-pos_tr)*0.01;
	fDepthLength = length(pend_tr-pos_tr);

	b = find_point_proxy(pos_tr, pend_tr, fStep, 5, 0.01, isec_tr);
	pos_tr = vec4(direction * pos);

	if (b)
	{
		fDepth = length(isec_tr-pos_tr)/fDepthLength;
		norm_tr = find_normal(isec_tr, 0.001);
		gl_FragColor = vec4((norm_tr.x+1.0)/2.0, (norm_tr.y+1.0)/2.0, (norm_tr.z+1.0)/2.0, fDepth);
	}
	else
	{
		gl_FragColor = vec4(0.0, 0.0, 0.0, 1.0);
	}
}
